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Dynamical reaction-diffusion processes and meta-population models are standard model- 
ing approaches for a wide variety of phenomena in which local quantities - such as den- 
sity, potential and particles - diffuse and interact according to the physical laws. Here, we 
study the behavior of two basic reaction-diffusion processes (B — > A and A + B —> 2B) 
defined on networks with heterogeneous topology and no limit on the nodes' occupation 
number. We investigate the effect of network topology on the basic properties of the sys- 
tem's phase diagram and find that the network heterogeneity sustains the reaction activ- 
ity even in the limit of a vanishing density of particles, eventually suppressing the critical 
point in density driven phase transitions, whereas phase transition and critical points, in- 
dependent of the particle density, are not altered by topological fluctuations. This work 
lays out a theoretical and computational microscopic framework for the study of a wide 
range of realistic meta-populations models and agent-based models that include the com- 
plex features of real world networks. 
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Reaction-diffusion (RD) processes are used to model phenomena as diverse as chemical re- 
actions, population evolution, epidemic spreading and many other spatially distributed systems 
in which local quantities obey physical reaction diffusion equations 12,3 . At the microscopic 
level, RD processes generally consist of particles (in many cases accounting for different kinds 
of "agents", information parcels, etc.) that diffuse in space and are subject to various reaction 
processes determined by the nature of the specific problem at hand. While fermionic RD mod- 
els assume exclusion principles that limit the number of particles on each node of the lattice, 
bosonic RD processes relax these constraints and allow each node of the lattice to be occu- 
pied by any number of particles. The classic example is provided by chemical reactions, in 
which different molecules or atoms diffuse in space and may react whenever in close contact. 
Another important instance for bosonic RD processes is found in meta-population epidemic 
models-^^&SJii. i n this case particles represent people moving between different locations, 
such as cities or urban areas. Individuals are divided into classes denoting their state with re- 
spect to the modeled disease — such as infected, susceptible, immune, etc. — and the reaction 
processes account for the possibility that individuals in the same location may get in contact 
and change their state according to the infection dynamics. 

The above modeling approaches are based on the spatial structure of the environment, trans- 
portation infrastructures, movement patterns, traffic networks, etc. The lack of accurate data on 
those features of the systems were usually reflected in the use of random graphs and regular 
lattices of different dimensionality as the substrate of the RD process. This corresponds to 
an implicit homogeneous assumption on the structure of the substrate, indeed used in many 
instances to solve the basic equations describing the RD process. In recent years, however, net- 
works which trace the activities and interactions of individuals, social patterns, transportation 
fluxes, and population movements on a local and global scale-^^iimj have been analyzed 
and found to exhibit complex features encoded in large scale heterogeneity, self-organization 
and other properties typical of complex systems 1 ^ 17 ' 18 ' 19 . In particular it has been found that a 
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wide range of societal and technological networks exhibit a very heterogeneous topology. The 
airport network among cities-!^, the commuting patterns in inter and intra-urban areas^ 1 ^, 
and several info-structures-^ are indeed characterized by networks whose nodes, representing 
the elements of the system, have a wildly varying degree, i.e. the number of connections to 
other elements. These topological fluctuations are mathematically encoded in a heavy-tailed 
degree distribution P(k), defined as the probability that any given node has degree k. They thus 
define highly heterogeneous substrates for the RD processes that cannot be accounted for in 
homogeneous or translationally invariant lattices. Analogously, models aimed at a description 
of spreading processes in spatially extended and societal systems are inevitably occurring in 
meta-population networks with connectivity patterns displaying very large fluctuations. Since 
connectivity fluctuations have been shown to have a large impact on the behavior of several 
percolation and fermionic systems 21122 ' 23 , the investigation of their role in the case of bosonic 
RD processes becomes a crucial issue for the understanding of a wide array of real world phe- 
nomena. 



Reaction-diffusion processes in complex networks. 

In order to investigate the effect of network heterogeneities on the phase diagram of meta- 
population models and chemical reaction processes, we consider a basic reaction scheme con- 
serving the number of particles that has been studied both in physics and mathematical epidemi- 
ology, namely the reaction-diffusion process identified by the following set of reactions 1 1 24 i 25 i 26 i 27 

B - A (1) 
B + A -> 2B. (2) 

From these reaction equations it is clear that the dynamics conserves the total number of par- 
ticles N = Na + Nb, where Ni is the number of particles i = A, B. This process can be 
naturally interpreted as a chemical reaction with an absorbing state phase transitio n 1 ' 25 ' 26 . The 
same reaction has been however used as a model problem in population dynamics in interaction 
with a polluting substance 27 , and it is analogous to the classic susceptible-infected-susceptible 
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(SIS) model for epidemic spreading 428 . In the process described by eqs. (QQ)-© the dynamics 
is exclusively due to B particles, that we identify as active particles, since A particles cannot 
generate spontaneously B particles. We consider the particles diffusing on a heterogeneous 
network with V nodes having a degree distribution P(k) characterized by the first and second 
moments (k) and (k 2 ), respectively. Reaction processes take place inside the network's nodes 
only, where each node i stores a number of A particles and hi of B particles (see Figure 1). 
The occupation numbers a; t and hi can assume any integer value, including a, t = hi = 0, that is, 
void nodes with no particles. For the sake of simplicity we assume that B particles diffuse with 
unitary time rate D B = 1 along one of the links departing from the node in which they are at 
a given time. This implies that at each time step a particle sitting on a node with degree k will 
jump into one of its nearest neighbor with probability 1/k. The results obtained in the following 
may be recovered for any diffusion rate D B ,at the expense of a more complicate mathematical 
treatment that will be reported elsewhere. In the case of A particles, we consider two different 
situations corresponding to a unitary (Da = 1) and a null (Da = 0) diffusion rate, respectively. 
While the first case is used in epidemic models which consider all individuals diffusing with 
the same rate, the case Da = is used in self-organized critical systems and specific absorbing 
phase transitions coupled with non-diffusive fields 28 i 25 i 26 . When Da = 0, the diffusion of A 
particles in the network occurs only through an effective process mediated by the reaction with 
B particles. Indeed any A particle may become a B particle following the reaction process and 
diffuse in the network until the reaction B — ► A occurs. This process is thus equivalent to an 
effective diffusion of A particles in the network. 

Before the diffusion process, the aj and hi particles stored in the same node react according 
to eqs. (OQ) and ©. In each node i the spontaneous process B — > A simply consists in turning 
each hi particle into an Oj particle with rate \i. We consider two general forms for the B + A — > 
2B process. In type I reaction we consider that each a, may react with all the hi particles in 
the same node, each reaction occurring with rate (3. In type II reaction we consider instead 
that each particle has a finite number of contacts with other particles. In this case the reaction 
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rate has to be rescaled by the total number of particles present in the node, i.e. /3/pi, where 
Pi = a,i + hi is the total number of particles in the node. This second process corresponds 
to what we usually observe in epidemic processes, where there is a population dependence of 
the reaction rate since individuals generally meet with a definite number of other individuals. 
In regular lattices and within the homogeneous mixing (mean-field) hypothesis, both type of 
processes exhibit a phase transition from an active phase (with an everlasting activity of B 
particles) to an absorbing phase (devoid of B particles), which in epidemic modeling correspond 
to the infected and healthy states, respectively. In the type I reaction the relevant parameter 
is represented by the average density of particles p = N/V and the transition occurs at the 
threshold value p = p c = p/P 25 ^ 2S . In the type II processes the transition point is found 
whenever (3/ p > 1. This second case is analogous to the classical epidemic threshold result and 
determines the lack or existence of an endemic state with a finite density of infected individuals 
(in this representation corresponding to the B particles) 34 . 

In order to take into account the topological fluctuations of the networks we have to explic- 
itly consider the presence of nodes with very different degree k. A convenient representation of 
the system is therefore provided by the quantities 

PA,k = — Y] CLi, p B ,k = — Y] h, (3) 

where Vk is the number of nodes with degree k and the sums run over all nodes % having de- 
gree ki equal to k. These two quantities express the average number of A and B particles in 
nodes with degree k. Analogously, p k = p A>k + p B ,k represents the average number of particles 
in nodes with degree k. The average density of A and B particles in the network is given by 
Pa = Y,k P(k)pA,k an d Pb = Y, k P(k)pB,k> respectively. Finally, by definition it follows that 
P = Pa + Pb- These quantities allow to express the RD process occurring on a heterogeneous 
network in terms of a set of rate equations describing the time evolution of the quantities pA,k(t) 
and p B) k{t) for each degree class k, as reported in the Materials and Methods section. The equa- 
tions depend on the reaction kernel T k that yields the number of B particles generated per unit 
time by the reaction processes taking place in nodes of a given degree class k. In uncorrelated 
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networks the resulting equations can be solved in the stationary limit, thus providing informa- 
tion on the phase diagram of the processes. 

Phase diagram and critical threshold. 

Let us first consider the type I reaction processes in uncorrelated networks. In this case we 
have that T k = pA,kPs,k obtaining in the stationary state 

P B,k = j^[(i-v)PB + pr], (4) 

where T = J2k P(k)T k . This equation readily states that the density of active particles is 
increasing in nodes with increasing degree k. This effect is mainly due to the diffusion process, 
that brings a large number of particles to well connected nodes, reflecting thus the impact of the 
network topological fluctuations on the particle density behavior (see the online supplementary 
information). In order to study the phase diagram we have to find the condition for which a 
solution p B 7^ of the set of equations for pA,k(t) and pB,k(t) is allowed. If D A = 0, i.e. for the 
case of non diffusing A particles, the density of A particles is independent of the node degree 
and is given by pA,k — Pa — p/P- hi view of the conservation of the number of particles this 
result readily implies that ps = p — ^ and therefore the presence of a phase transition from an 
absorbing phase to an active state at a critical value of the total density of particles p c — <|. A 
very different picture is obtained when A particles are also allowed to diffuse. In this case, for 
D A = 1, the stationary density of A particles is given by 

k 

pA,k = j^(pA + PPb - I3T) . (5) 

The system of equations can be solved by imposing a self-consistent condition for the quantity 
T (the details of the calculation are reported in the online supplementary information) and the 
non-trivial solution p B > is obtained only if the total density of particles satisfies the condition 
p > Pc, with 
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This result implies that, if A particles can diffuse, topological fluctuations affect the critical 
value of the transition. Networks characterized by heterogeneous connectivity patterns display 
large degree fluctuations so that (k 2 ) ^> (k) 2 . In the infinite size limit V — > oo we have (A; 2 ) — > 
oo and thus eq. © yields a critical value p c = 0, showing that the topological fluctuations of 
the network suppress the phase transition in the thermodynamic limit. This is a very relevant 
result that, analogously to those concerning percolation- 2 ^ 2 , and standard epidemic processes 23 , 
indicates that physical and dynamical processes taking place on scale-free and heavy-tailed 
networks behave very differently with respect to the same processes occurring on homogeneous 
networks. 

In Figure 2 we provide support to this theoretical picture by reporting the results obtained 
from Monte Carlo simulations of RD processes of type I on uncorrected networks with given 
scale-free degree distribution P(k) ~ /c~ 7 . The simulations use a single particle modeling 
strategy in which each individual particle is tracked in time. The system evolves following a 
stochastic microscopic dynamics and at each time step it is possible to record average quanti- 
ties, such as e.g. the density of active particles ps{t). In addition, given the stochastic nature 
of the dynamics, the experiment can be repeated with different realizations of the noise, differ- 
ent underlying graphs, and different initial conditions. This approach is equivalent to the real 
evolution of the RD process in the generated networks and can be used to validate the theo- 
retical results obtained in the analytical approach. The top panel of Figure 2 shows the phase 
transitions observed in the two cases, whether A particles diffuse or not. If Da = the process 
undergoes a phase transition at p c = p,/ (3 = 2, regardless of the difference in the level of hetero- 
geneity as provided by different values of the power-law exponent 7 of the degree distribution 
P(k) at fixed network size (here V = 10 4 , 7 = 3 and 7 = 2.5) or by different sizes V at fixed 
7 (7 = 2.5, V = 10 4 and V = 10 5 ). In case A particles diffuse, instead, the transition occurs at 
critical values p c < p//3, with p c — > for decreasing ratios (k) 2 / (k 2 ) as observed for increasing 
sizes V at fixed 7 (curves for 7 = 2.5 with V from 10 3 to 10 5 nodes are shown), in agreement 
with the analytical result of eq. ©. Bottom panels show the difference in the behavior of pA,k 



7 



Colizza et al. 



as a function of the degree k: a flat spectrum is obtained when D A = and a linear dependence 
in k when Da = 1 • 

A different scenario emerges when considering type II processes. In this case the reaction 
kernel is Tk = PA,kPB,k/ Pk\ i- e - m eacn node A particles will participate in a reaction event 
with a rate proportional to the relative density of B particles. While the set of equations for 
pB,k has the same form of eq. ©, the stationary condition for pA,k yields solutions that depend 
on k for both diffusive and non-diffusive A particles. In particular, in both cases we have that 
Pb > if the condition (3/ p > 1 is satisfied (see the online supplementary information). This 
result recovers the usual threshold condition which depends only on the reaction rates and is not 
affected by changes in the total density of particles p. Also for type II processes we performed 
extensive Monte-Carlo simulations considering uncorrelated scale-free networks with a heavy 
tailed degree distribution. Figure 3 reports the results obtained in the two cases Da = and 
Da = 1, with different underlying network topologies. Changes in the number of nodes and 
in the exponent 7 assumed for the degree distribution do not affect the phase transition. The 
critical value depends exclusively on the process rates, despite the observed linear behavior of 
PA,k and p B ,k bears memory of the heterogeneity of the underlying network. 

Discussion and comparison with realistic models. 

The different phase diagrams obtained in type I and II processes can be understood qualita- 
tively in terms of the following argument. In type I processes, whatever the parameters (3 and 
p, there exists a value of the local density large enough to keep the system in an active state 
by sustaining the creation of B particles in the right amount. Large topological fluctuations 
imply the existence of high degree nodes with a high density of particles and therefore a high 
number of generated B particles. This implies that in the thermodynamic limit there is always 
a node (with a virtually infinite degree) with enough particles to keep alive the process even for 
a vanishing average density of particles, leading to the suppression of the phase transition. The 
crucial ingredient of this mechanism is given by the diffusion process that allows high degree 
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nodes to have a number of particles that is proportional to their degree. This is confirmed by 
the case in which A particles do not diffuse. In this case high degree nodes do not accumulate 
enough particles and the usual threshold effect is recovered. In type II processes, on the other 
hand, the reaction activity in each node is rescaled by the local density pi and it is therefore the 
same in all nodes, regardless of the local population. In this case, the generation of B particles 
is homogeneous across nodes of different degrees and therefore the presence of an active state 
depends only on the balance between the reaction rates (3 and //. 

These results let emerge a basic framework for the microscopic (mechanistic in the epi- 
demic terminology) description of meta-population epidemic models. The type I and II pro- 
cesses correspond to the two limits of transmissibility independent or inversely proportional to 
the population size, respectively. In addition, realistic meta-population models have heteroge- 
neous diffusion probabilities due to the traveling pattern and fixed population sizes according 
to data. Despite these extra complications, the basic reaction-diffusion framework studied here 
provides a simple qualitative picture of the realistic models. In Figure 4 we report the two types 
of processes studied here and compare them with the results from a realistic compartmental 
SIS meta-population model considering 500 urban areas in the Unites States and including the 
actual data of the air traveling flows among those urban area s 33134 . The network is defined by 
nodes representing each urban area together with its population and edges representing air travel 
fluxes along which individuals diffuse, coupling the epidemic spreading in different urban areas 
(see Ref. 3 ^^ for a detailed definition of the model). The type I process is compared with a 
model whose transmissibility is independent of the population size and the type II process is 
compared with the usual epidemic spreading with trasmissibility scaling proportionally to the 
population size (see Ref. 5 ). Figure 4 shows in the four cases the reaction activity occurring on 
each network node, i.e. the creation of B particles in the reaction-diffusion process and newly 
infected individuals in the realistic epidemic model normalized to the local population. Increas- 
ing values of the reaction activity correspond to colors ranging from yellow (low activity) to 
red (high activity). In type I processes the reaction activity is linearly increasing with the pop- 
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ulation of the nodes, thus showing high activity (red) concentrated in largely populated nodes 
(represented with a larger size). The homogeneity in the generation of B particles in type II 
processes is evident: all nodes display the same color and thus experience the same level of ac- 
tivity, regardless of the local density. Strikingly, despite the various complications and elements 
of realism introduced in the data-driven meta-population model, its qualitative behavior is in 
very good agreement with the results obtained for the microscopic reaction-diffusion processes 
in both transmissibility limits. 

In summary, the microscopic RD framework introduced here is able to provide a general 
theoretical understanding of the behavior of more realistic meta-population epidemic mod- 
els. Furthermore, the presented approach can be extended to include the various sources of 
heterogeneity — such as degree correlations-^^, heterogeneous diffusion probabilities and their 
non-linear relations with the connectivity pattern — needed in order to provide a detailed analy- 
sis of realistic processes. 

Materials and Methods 

Reaction- diffusion equation. In order to take fully into account degree fluctuations in an 
analytical description of the RD processes, we have to relax the homogeneity assumption and 
allow for degree fluctuations by introducing the relative densities pB,k(t), PA,k(t), and pk(t). 
The dynamical reaction rate equations for B particles in any given degree class can thus be 
written as 

d tP B,k = ~PB,k + kY,P(k'\k)^[(l - p) PB ,y + /?IV] , (7) 
k> K 

where P(k'\k) represents the conditional probability that a vertex of degree k is connected to 
a vertex of degree k '— . The various terms of the equations are obtained by considering that at 
each time step the particles present on a node of degree k first react and then diffuse away from 
the node with unitary diffusion rate accounted by the term —ps,k- The positive contribution for 
the particle density is obtained by summing the contribution of all particles diffusing in nodes 
of degree k from their neighbors of any degree k' , including the new particles generated by the 
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reaction term T k . In the case of uncorrelated networks the conditional probability P{k'\k) that 
any given edge points to a vertex with k' edges is independent of k and equal to k'P(k') / {k)^2., 
so that the reaction rate equations read as 

k 

d t p B ,k = ~pB,k + -rjT [(1 - P)Pb + PT] , (8) 

where pB = Y,k P(k)pB,k and T = Y,k P{k)Y k . In the case for pA,k(t) we have to distinguish 
whether if A particles diffuse or not. If Da = 1, we obtain a set of equations analogous to those 
for p B ,k that read as 

k 

d t pA,k = ~PA,k + jj^(pa + PPb - (3T) , (9) 

where pa = J2k P{k)pA,k- In the case of non-diffusive A particles (Da = 0) the equations 
reduce to: 

d t PA,k = PPB,k - (3Tk- (10) 

The phase diagram for the various cases and the conditions for p B > are obtained by imposing 
the stationary state defined by d t pA,k = and d t ps,k = 0, with the additional constraint that 
P = Pa + Pb, i- e - the number of particles is conserved. We are therefore led to a simple set of 
algebraic equations whose explicit solution is reported in the online supplementary information. 

Monte-Carlo simulations. The uncorrelated networks considered have been generated with 
the uncorrelated configuration model 31 , based on the Molloy-Reed 32 algorithm with an ad- 
ditional constraint on the possible maximum value of the degree in order to avoid inherent 
structural correlations. The algorithm is defined as follows. Each node i is assigned a degree 
ki obtained from a given degree sequence P(k) (in our case P(k) ~ fc~ 7 with 7 = 3 and 
7 = 2.5) subject to the restriction ki < V 1 / 2 . Links are then drawn to randomly connect pairs 
of nodes, respecting their degree and avoiding self-loops and multiple edges. Sizes of V = 10 3 , 
V = 10 4 and V = 10 5 nodes have been considered. Initial conditions are generated by ran- 
domly placing Vpa(0) particles A and Vps(0) particles B, corresponding to a particle density 
P = Pa(0) + Pb(0). The results are independent of the particular initial ratio pa(0)/ Pb(0), 
apart from very early time transients. The dynamics proceeds in parallel and considers discrete 
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time steps representing the unitary time scale r of the process. The reaction and diffusion rates 
are therefore converted into probabilities and at each time step, the system is updated accord- 
ing to the following rules, a) Reaction processes: i) On each lattice site, each B particle is 
turned into an A particle with probability pr. ii) At the same time, each A particle becomes a 
B particle with probability determined by the type of reaction process, b) After all nodes have 
been updated for the reaction, we perform the diffusion: on each lattice site, each B particle 
moves into a randomly chosen nearest neighbor site; the same process occurs for A particles if 
D A = 1. The simulation details of the reaction process represented by eq. © depend on the 
kernel considered. In type I processes each A particle in a given node i becomes a B particle 
with probability 1 — (1 — f3r) h \ where bi is the total number of B particles in that node. This 
corresponds to the average probability for an A particle of being involved in the reaction © 
with any of the B particles present on the same site. In type II processes the reaction process 
is simulated by turning each A particle into a B particle with probability 1 — (1 — j;r) bl where 
Pi is the total number of particles in the node i. This term accounts for the average probability 
that an A particle will get in contact with a B particle present in the node, given that the pos- 
sible number of contacts is rescaled by the population p, of the node. The term (3/ therefore 
represents the normalized transmission rate of the process. 

Correspondence and requests for material should be addressed to V.C. (vcolizza@indiana.edu) 
or A.V. (alexv@indiana.edu). 
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Figure captions 

Figure 1. Bosonic reaction-diffusion systems in heterogeneous networks. Schematic repre- 
sentation of RD processes in heterogeneous complex networks when the multiple occupancy of 
nodes is allowed. Particles A and B can diffuse in the network and, inside each node, undergo 
the reaction processes described by eqs. CD-©- Each node i stores pi = + hi particles, where 
the occupation numbers ai and b{ can assume any integer value, including zero. 

Figure 2. Phase diagram and stationary densities for type I processes. Top panel: Phase 
transitions in type I processes for diffusing and non-diffusing A particles. If Da = the 
transition occurs at the critical value of the density p c = p/ (3 = 2, regardless of the topol- 
ogy of the underlying network. Results for uncorrected scale-free networks having degree 
distribution P(k) ~ fc -7 with 7 = 2.5 and 7 = 3.0 and different sizes V show the same 
behavior, small differences in the value of ps being due to finite- size effects. If Da = 1, 
the critical point is strongly affected by the topological fluctuations of the network. Here 
we show results for 7 = 2.5 and sizes of the network V = 10 3 , 10 4 , 10 5 corresponding to 
(k) 2 /(k 2 ) = 0.52, 0.32, 0.19, respectively. With increasing sizes, degree fluctuations become 
larger and the transition is observed at smaller values of p c . Bottom panels: Stationary densities 
pA,k and pB,k as functions of the degree k. If Da = 0, the average density of A particles inside 
nodes of degree k is constant, while the behavior shown by B particles is linear in k. If Da = 1, 
both densities are linear in k. 

Figure 3. Phase diagram and stationary densities for type II processes. Top panel: Phase 
transitions in type II processes for diffusing and non-diffusing A particles. Regardless of topo- 
logical fluctuations in the underlying network and of the probability of diffusion Da, the tran- 
sition occurs at the critical point (3/ p = 1, depending only on the reaction rates. Here we show 
results for networks of size V = 10 4 with particle density p = 20, power-law exponents 7 = 2.5 
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and 7 = 3.0. Differences in the values of the stationary density p B are due to finite- size effects. 
Bottom panels: Stationary densities pA.k and ps,k as functions of the degree k. In both cases, 
Da = and Da = 1, linear dependencies in k are obtained. 

Figure 4. Reaction activity in type I and type II processes: microscopic model and real- 
world examples. The first row of panels refers to the microscopic RD model, as described in the 
text, while the second row reports the analysis of the spread of an airline-carried disease in the 
US with a data-driven meta-population model. Both models consider the actual topology of the 
US air transportation network as obtained by considering the 500 airports with largest traffic 33 ; 
nodes population is obtained from census data 34 . In addition, the realistic meta-population 
model also considers the actual traffic of passengers on each connection between airports. The 
networks are mapped on a globe for sake of visualization. Each node is represented with a 
size linearly dependent on its population and a color illustrating the level of reaction activity 
inside the node, ranging from to the max value. Whereas type I processes experience a level 
of activity proportional to nodes population - corresponding to red color in largely populated 
nodes and yellow in small population nodes - the reaction activity is homogeneously distributed 
among the nodes of the network when type II processes are considered. 
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Figure 1. 
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Figure 3. 
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